Figure 1 plots

Load libraries

Loading required package: Matrix
Loading required package: igraph
package ‘igraph’ was built under R version 3.5.3
Attaching package: ‘igraph’

The following objects are masked from ‘package:stats’:

    decompose, spectrum

The following object is masked from ‘package:base’:

    union


Attaching package: ‘dplyr’

The following objects are masked from ‘package:igraph’:

    as_data_frame, groups, union

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


Attaching package: ‘conos’

The following objects are masked _by_ ‘.GlobalEnv’:

    plotClusterBarplots, plotComponentVariance, plotDEheatmap

Loading required package: ggplot2
Learn more about the underlying theory at https://ggplot2-book.org/

Attaching package: ‘cowplot’

The following object is masked from ‘package:ggplot2’:

    ggsave

Custom clustering

#ncdcon$findCommunities(method=rleiden.community,r=c(1,0.5))
#ncdcon$findCommunities(method=leiden.community,r=1.2)

Load NB-only alignment

ncdcon <- readRDS("~pkharchenko/m/ninib/NB/ncdcon.rds")

doubletScores <- unlist(readRDS("~pkharchenko/m/ninib/NB/doubletScores.rds"))

Initial annotation

ninib_bannot <- readRDS("~pkharchenko/m/ninib/NB/ninib_bannot.rds")
annot <- readRDS("/d0-mendel/home/meisl/Workplace/neuroblastoma/cell.annotation.Jan2020.rds")$cellano

Clean up annotation

annot <- readRDS("/d0-mendel/home/meisl/Workplace/neuroblastoma/cell.annotation.Jan2020.rds")$cellano
annot <- setNames(as.character(annot),names(annot))
annot[annot=='unknown'] <- 'Th' # Ca2+ mediated apoptosis, involivng LCK, CD3D, RAC2, etc. Lumping in with Th 
annot[annot=='Plasmacytoid'] <- 'pDC'
annot[grep('Cytotoxic',annot)] <- 'Tcyto'
annot[grep('T helper',annot)] <- 'Th'
annot[grep('NK',annot)] <- 'NK'
annot[grep('Bcell',annot)] <- 'B'
annot[grep('PlasmaCell',annot)] <- 'Plasma'
annot[grep('Mast',annot)] <- 'Mast'
annot[grep("Bridge",annot)] <- 'SCP-like'
dannot <- as.factor(annot)

annot[grep("SOX11|Stress|euronal|Prolif",annot)] <- 'Noradrenergic'
annot[grep("Bridge",annot)] <- 'SCP-like'
annot <- as.factor(annot);
# renamings
levels(annot) <- gsub("Noradrenergic","Adrenergic",levels(annot))

set up pallets

set.seed(3)
annot.pal <- setNames(sample(rainbow(length(levels(annot)),v=0.95,s=0.7)),levels(annot));
annot.palf <- function(n) return(annot.pal)
dannot.pal <- setNames(sample(rainbow(length(levels(dannot)),v=0.85,s=0.9)),levels(dannot));
x <- which(levels(dannot) %in% levels(annot)); dannot.pal[x] <- annot.pal[match(levels(dannot)[x],levels(annot))]
dannot.palf <- function(n) return(dannot.pal)

# for samples
sample.pal <- setNames(sample(rainbow(length(ncdcon$samples),v=0.7,s=0.9)),names(ncdcon$samples))
sample.pal <- sample.pal[c(sort(names(sample.pal))[-1],'Adr')]
sample.palf <- function(n) return(sample.pal[1:n])
sample.palf <- function(n) { browser(); return(sample.pal[1:n])}

p1 <- ncdcon$plotGraph(alpha=0.03,size=0.5,groups=ncdcon$getDatasetPerCell(),plot.na=F,mark.groups=F,palette=sample.pal)
p2 <- ncdcon$plotGraph(alpha=0.03,size=0.5,groups=annot,plot.na=F,mark.groups=T,palette=annot.palf)
#p3 <- ncdcon$plotGraph(alpha=0.03,size=0.5,plot.na=F,mark.groups=T,clustering='leiden')
p3 <- ncdcon$plotGraph(alpha=0.03,size=0.5,groups=dannot,plot.na=F,mark.groups=T,palette=dannot.palf)
plot_grid(plotlist=list(p1,p2,p3),nrow=1)

write out p2 app

source("~/m/pavan/DLI/conp2.r")
cell.subset <- unique(c(sample(rownames(ncdcon$embedding),2e3),names(annot)[annot=='SCP-like']))
ncdcon.p2 <- p2app4conos(ncdcon,file='ncdcon.bin',metadata = list(coarse=annot,detailed=dannot,cnv=cnv.sc),cell.subset = cell.subset)

Sample panel:

pp <- ncdcon$plotPanel(groups=annot,plot.na=F,alpha=0.1,size=0.2,use.common.embedding = T, ncol=2, mark.groups=F,palette=annot.pal,raster=T, raster.width=1,raster.height=1,title.size=4)
#pdf('panel.pdf',height=8,width=2); print(pp); dev.off();
pp

ncdcon$plotGraph(gene='TOP2A',alpha=0.5,size=0.1)

ncdcon$plotGraph(gene='TOP2A',alpha=0.5,size=0.1)

source("~/m/pavan/DLI/conp2.r")
ncdcon.p2 <- p2app4conos(ncdcon,file='ncdcon.bin',max.cells=1e4,metadata = list(annot=as.factor(annot),dannot=as.factor(dannot)))

Marker genes:

nfac <- annot;

nfac.de <- ncdcon$getDifferentialGenes(groups=nfac,n.cores=30,append.auc=TRUE,z.threshold=0,upregulated.only=T)

source("~/m/p2/conos/R/plot.R")
pp <- plotDEheatmap(ncdcon,nfac,nfac.de,n.genes.per.cluster = 20 ,show.gene.clusters=T,column.metadata=list(samples=ncdcon$getDatasetPerCell()),order.clusters = T,use_raster=F, column.metadata.colors = list(clusters=annot.pal,samples=sample.pal))
pdf(file='nfac.heatmap2.pdf',width=10,height=30); print(pp); dev.off();

A small version of the hetmap, for the main figure

source("~/m/p2/conos/R/plot.R")
genes <- c("CD79B",'MS4A1','PLVAP','EPCAM','SPINK2','IL7R','LYZ','C1QB','CLEC10A','COL1A1','LUM','S100A9','IL1B','MYL9','GNLY','RTN1','MDK','NNAT','IRF4','MYH11','FCRL5','PLP1','GZMA','RORA','TIGIT')
pp <- plotDEheatmap(ncdcon,nfac,nfac.de,n.genes.per.cluster = 20 ,show.gene.clusters=T,column.metadata=list(samples=ncdcon$getDatasetPerCell()), column.metadata.colors = list(clusters=annot.pal,samples=sample.pal), order.clusters = T, additional.genes = genes, labeled.gene.subset = genes, use_raster=F,min.auc = 0.79)
pp

pdf(file='markers.small.pdf',width=6,height=6); print(pp); dev.off();

Make a very simplified version

sannot <- setNames(as.character(annot),names(annot));
sannot[grep("^T|ILC3|NK",sannot)] <- 'T cells'
sannot[grep("^B$|Plasma|pDC",sannot)] <- 'B cells'
sannot[grep("Macrop|Mono|mDC|Mast",sannot)] <- 'Myeloid'
sannot <- as.factor(sannot)

sannot.pal <- setNames(annot.pal[match(levels(sannot),names(annot.pal))],levels(sannot))
sannot.pal['T cells'] <- annot.pal['Tcyto']
sannot.pal['B cells'] <- annot.pal['B']
sannot.pal['Myeloid'] <- annot.pal['Macrophages']

sannot.de[['Noradrenergic']] %>% arrange(-AUC) %>% head(30)

sannot.de[['Adrenergic']] %>% arrange(-AUC) %>% head(30)

source("~/m/p2/conos/R/plot.R")
pp <- plotDEheatmap(ncdcon,sannot,sannot.de,n.genes.per.cluster = 20 ,show.gene.clusters=T,column.metadata=list(samples=ncdcon$getDatasetPerCell()),order.clusters = T, column.metadata.colors = list(clusters=sannot.pal,samples=sample))
pdf(file='sannot.heatmap.pdf',width=10,height=30); print(pp); dev.off();

A small version of the hetmap, for the main figure

source("~/m/p2/conos/R/plot.R")
genes <- c("CD79B",'MS4A1','PLVAP','EPCAM','SPINK2','IL7R','LYZ','C1QB','CLEC10A','COL1A1','LUM','S100A9','IL1B','MYL9','GNLY','RTN1','MDK','NNAT','IRF4','MYH11','FCRL5','PLP1','GZMA','RORA','TIGIT')
genes <- c("CD3D",'CD7','CD79A','HLA-DRA','CD74','LYZ','STMN2','MDK','MYL9','MYH11','PECAM1','LUM','PLP1')
tf <- sannot[unlist(tapply(1:length(sannot),sannot,function(x) { if(length(x)>2e2) x <- sample(x,2e2); x}))]
pp <- plotDEheatmap(ncdcon,tf,sannot.de,n.genes.per.cluster = 20 ,show.gene.clusters=T,column.metadata=list(samples=ncdcon$getDatasetPerCell()), column.metadata.colors = list(clusters=sannot.pal,samples=sample.pal), order.clusters = T, additional.genes = genes, labeled.gene.subset = genes, min.auc = 0.75,use_raster = T,raster_device = "CairoPNG",averaging.window=5)
pp

pdf(file='sannot.select.pdf',width=4,height=6); print(pp); dev.off();

source("~/m/p2/conos/R/plot.R")
pdf(file='sannot.select.pdf',width=5,height=6); 
plotDEheatmap(ncdcon,sannot,sannot.de,n.genes.per.cluster = 20 ,show.gene.clusters=T,column.metadata=list(samples=ncdcon$getDatasetPerCell()), column.metadata.colors = list(clusters=sannot.pal,samples=sample.pal), order.clusters = T, additional.genes = genes, labeled.gene.subset = genes, min.auc = 0.75)
dev.off();

pdf(file='sannot.small.pdf',width=6,height=6); print(pp); dev.off();

null device 
          1 

pdf(file='sannot.small.pdf',width=6,height=6); print(pp); dev.off();

Draw just the tumor cell annotations:

source("~/m/p2/conos/R/plot.R")
tf <- droplevels(sannot[grep("Adr|SCP|Mesench",sannot)]); tf <-factor(tf,levels=levels(tf)[c(2,3,1)])
pp <- plotDEheatmap(ncdcon,tf,sannot.de,n.genes.per.cluster = 10, expression.quantile=0.95, show.gene.clusters=T,column.metadata=list(samples=ncdcon$getDatasetPerCell()),order.clusters = F, column.metadata.colors = list(clusters=sannot.pal,samples=sample.pal),additional.genes = c("SOX10","PRRX1", "LEPR", "PDGFRA","TH","ISL1","NRG1","ERBB3"),exclude.genes=c("CRYAB"),min.auc=0.75,row.label.font.size = 10,split=T, split.gap=1,column_title = NULL,use_raster = T,raster_device = "CairoPNG",averaging.window=5)
pp

pdf(file='NB.markers.small.pdf',width=4,height=6); print(pp); dev.off();

null device 
          1 

ncdcon$plotGraph(alpha=0.3,size=0.5,plot.na=F,mark.groups=T,gene='RORA')

Plot CNV scores

ncdcon$plotPanel(gene='PLP1',alpha=1,size=0.3, use.common.embedding = T,ncol=2,gradient.range.quantile=0.8)

See interaction between PLP1 expression in the ‘bridge’ cell population and the CNV scores across different samples

#exp <- pmin(conos:::getGeneExpression(ncdcon,'PLP1'),conos:::getGeneExpression(ncdcon,'S100B'))

exp <- conos:::getGeneExpression(ncdcon,'PLP1')
samplef <- ncdcon$getDatasetPerCell()

pl <- lapply(names(cutoff.list),function(san) {
  cn <- intersect(names(annot)[annot=='SCP-like'],names(samplef)[samplef==san])
  df <- data.frame(cnv=cnv.sc[match(cn,names(cnv.sc))],exp=exp[match(cn,names(exp))],row.names=cn)
  #ggplot(df,aes(exp,cnv)) + geom_point(col=adjustcolor(1,alpha=0.3)) + ggtitle(san) +geom_hline(yintercept = 0, linetype="dashed", color = "red")
  #ggplot(df,aes(exp,cnv)) + geom_point(col=adjustcolor(pagoda2:::val2col(doubletScores[cn],zlim=c(0,0.5),gradientPalette=colorRampPalette(c('black','red'))(1024)),alpha=0.3)) + ggtitle(san) +geom_hline(yintercept = 0, linetype="dashed", color = "red")
  ggplot(df,aes(exp,cnv)) + geom_point(col=adjustcolor(ifelse(doubletScores[match(cn,names(doubletScores))]>0.2,'red','black'),alpha=0.3)) + ggtitle(san) +geom_hline(yintercept = 0, linetype="dashed", color = "red")
})

plot_grid(plotlist=pl)

Just the scores

cnv.sc <- readRDS("/d0-mendel/home/meisl/Workplace/neuroblastoma/AmpScore.rds")

Stringent scores:

load("/d0-mendel/home/meisl/Workplace/neuroblastoma/Figures/S1/CNV/FDR5.RData")
names(cutoff.list) <- gsub("\\..*","",names(cutoff.list))
cutoff.list <- setNames(unlist(cutoff.list),names(cutoff.list))

valid.samples <- c('NB12','NB26')
#valid.samples <- unique(gsub("_.*","",names(cnv.sc)));
cnv.sc <- cnv.sc[ncdcon$getDatasetPerCell()[names(cnv.sc)] %in% valid.samples]

Detailed info allScoreA: scaled gene set average expression. cutoff.list: cutoff for each patient CNV.cells: amplified cell name

load("/d0-mendel/home/meisl/Workplace/neuroblastoma/Figures/CNV/CNV.0211.RData")
names(cutoff.list) <- gsub("\\..*","",names(cutoff.list))
cutoff.list <- setNames(unlist(cutoff.list),names(cutoff.list))

Calculate a scaled score

cell.cnv.thresholds <- setNames(cutoff.list[gsub("_.*","",names(allScoreA))],names(allScoreA))
cnv.sc <- allScoreA-cell.cnv.thresholds

Limit to the samples where we see tumor cells and CNVs

valid.samples <- c('NB01','NB09','NB12','NB13','NB15','NB18','NB19','NB20','NB21','NB24','NB26')
cnv.sc <- cnv.sc[ncdcon$getDatasetPerCell()[names(cnv.sc)] %in% valid.samples]

hist(cnv.sc,col='wheat')

emb <-ncdcon$embedding
ncdcon$embedding <- emb[names(cnv.sc)[order(cnv.sc)],]
col <- pagoda2:::val2col(cnv.sc,gradient.range.quantile = 0.9,gradientPalette = colorRampPalette(c(rep('grey85',3),'red'))(1024))
p3 <- ncdcon$plotGraph(alpha=0.1,size=0.5,colors=col,plot.na=F)
ncdcon$embedding <- emb;
p3

Threshold-based, with overplotting

col <- cnv.sc>0.7;
#col <- cnv.sc>0.8e-3; 
emb <-ncdcon$embedding
ncdcon$embedding <- emb[names(col)[order(col)],]
p3 <- ncdcon$plotGraph(alpha=0.05,size=0.6,groups=col,mark.groups=F,raster=T,palette=c("TRUE"='red',"FALSE"='gray90'),plot.na=F,raster.height=4,raster.width=4)
ncdcon$embedding <- emb;
#pdf(file='cnv_overview.pdf',width=4,height=4); print(p3); dev.off();
p3

Combined figure

p1 <- ncdcon$plotGraph(alpha=0.03,size=0.5,groups=ncdcon$getDatasetPerCell(),plot.na=F,mark.groups=F,palette=sample.pal)
p2 <- ncdcon$plotGraph(alpha=0.03,size=0.5,groups=annot,plot.na=F,mark.groups=T,palette=annot.palf)
#p3 <- ncdcon$plotGraph(alpha=0.03,size=0.5,plot.na=F,mark.groups=T,clustering='leiden')
#p3 <- ncdcon$plotGraph(alpha=0.03,size=0.5,colors=pagoda2:::val2col(cnv.sc,gradient.range.quantile = 0.9),plot.na=F)
pl <- list(p1,p2,p3)
plot_grid(plotlist=pl,nrow=1)

require(gridExtra)

Loading required package: gridExtra

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

raster <- T; size=0.5; alpha=0.03;
p1 <- ncdcon$plotGraph(alpha=alpha,size=size,raster=raster,groups=ncdcon$getDatasetPerCell(),plot.na=F,mark.groups=F,palette=sample.pal,raster.height=5,raster.width=5)
p2 <- ncdcon$plotGraph(alpha=alpha,size=size,raster=raster,groups=annot,plot.na=F,mark.groups=T,palette=annot.palf,font.size=c(3,5),raster.height=5,raster.width=5)
#p3 <- ncdcon$plotGraph(alpha=alpha,size=size,raster=raster,colors=col,plot.na=F)
pl <- list(p1,p2,p3)
pdf(file='overview.pdf',width=4,height=6); 
grid.arrange(p2,p1,p3,layout_matrix=rbind(c(1,1),c(2,3)),heights=c(2,1))
dev.off();

null device 
          1 

invisible(lapply(1:length(pl),function(i) { pdf(paste('panel',i,'pdf',sep='.'),width=4,height=4); print(pl[[i]]); dev.off()}))

Looking at cells with SCP-like pluripotency genes

ascl1e <- conos:::getGeneExpression(ncdcon,'ASCL1')
phox2be <- conos:::getGeneExpression(ncdcon,'PHOX2B')
prrx2e <- conos:::getGeneExpression(ncdcon,'PRRX2')
prrx1e <- conos:::getGeneExpression(ncdcon,'PRRX1')
sox10e <- conos:::getGeneExpression(ncdcon,'SOX10')
foxd3e <- conos:::getGeneExpression(ncdcon,'FOXD3')

dc <- sox10e>0 & phox2be>0

#dc <- sox10e>0 & ascl1e>0
dc <- sox10e>0 & phox2be>0
dc <- foxd3e>0 & phox2be>0
dc <- sox10e>0 & prrx2e>0
dc <- sox10e>0 & prrx1e>0
dc <- phox2be>0 & prrx1e>0
dc <- phox2be>0.3 & prrx1e>0.3 & sox10e>0
table(dc)
names(dc)[dc]
em <- ncdcon$samples$NB26$embeddings$PCA[[2]]
plot(em[,1],em[,2],pch=19,cex=0.2,col=adjustcolor(1,alpha=0.01))
vc <- names(dc)[dc]
points(em[rownames(em) %in% vc,1],em[rownames(em) %in% vc,2],col=2)

Look at scatter plots of key fork markers

gns <- c('SOX10','ASCL1','PHOX2B','PRRX2','PRRX1','FOXD3','POSTN')
gnse <- lapply(pagoda2:::sn(gns),function(g) conos:::getGeneExpression(ncdcon,g))

g1 <- 'SOX10'; g2 <-'PHOX2B'
g1 <- 'POSTN'; g2 <-'PHOX2B'
#g1 <- 'SOX10'; g2 <-'PRRX2' # 1 cell
#g1 <- 'SOX10'; g2 <-'PRRX1' # clear
#g1 <- 'SOX10'; g2 <-'ASCL1' # none
#g1 <- 'FOXD3'; g2 <- 'PHOX2B'
#g1 <- 'PRRX1'; g2 <- 'PHOX2B'

dc <- gnse[[g1]]>0 & gnse[[g2]]>0
vc <- names(dc)[dc]; vc <- vc[grepl('^NB',vc)]
df <- data.frame(gnse[[g1]][vc],gnse[[g2]][vc]); colnames(df) <- c(g1,g2)
df <- cbind(df,data.frame(type=annot[vc],cell=vc,doublet=doubletScores[vc]>0.2,inemb=vc %in% rownames(ncdcon$embedding)))
p1 <- ggplot(df,aes_(x=as.name(g1),y=as.name(g2),colour=quote(type))) + geom_point(size=2) + scale_color_manual(values=annot.pal)
p2 <- ggplot(df,aes_(x=as.name(g1),y=as.name(g2),colour=quote(doublet),shape=quote(inemb))) + geom_point() + scale_color_manual(values=c("FALSE"="gray30", "TRUE"="red"));
cowplot::plot_grid(plotlist=list(p1,p2))

Make a panel of pairs

pl <- list(c('SOX10','PHOX2B'),c('SOX10','PRRX1'),c('SOX10','PRRX2'),c('PRRX1','PHOX2B'),c('POSTN','PHOX2B'))
pp <- lapply(pl,function(x) {
  g1 <- x[1]; g2 <- x[2]
  dc <- gnse[[g1]]>0 & gnse[[g2]]>0
  vc <- names(dc)[dc]; vc <- vc[grepl('^NB',vc)]
  df <- data.frame(gnse[[g1]][vc],gnse[[g2]][vc]); colnames(df) <- c(g1,g2)
  df <- cbind(df,data.frame(type=annot[vc],cell=vc,doublet=doubletScores[vc]>0.2,inemb=vc %in% rownames(ncdcon$embedding)))
  p1 <- ggplot(df,aes_(x=as.name(g1),y=as.name(g2),colour=quote(type))) + geom_point(size=2) + scale_color_manual(values=annot.pal)
})
cowplot::plot_grid(plotlist=pp,ncol=2)

vc <- names(dc)[dc]; vc <- vc[grepl('^NB',vc)]
df <- data.frame(phox2b=phox2be[vc],prrx1=prrx1e[vc],type=annot[vc],cell=vc)
ggplot(df,aes(x=phox2b,y=prrx1,colour=type)) + geom_point()

dc <- phox2be>0 & prrx1e>0
dc <- phox2be>0 & prrx1e>0
dc <- phox2be>0.3 & prrx1e>0.3 & sox10e>0
vc <- names(dc)[dc]; vc <- vc[grepl('^NB',vc)]

p1 <- ncdcon$plotGraph(alpha=0.03,size=0.5,groups=annot,plot.na=F,mark.groups=T)
p2 <- ncdcon$plotGraph(alpha=0.03,size=0.5,gene='SOX10',plot.na=F,mark.groups=T)
p3 <- ncdcon$plotGraph(alpha=0.2,size=3,groups=setNames(rep(1,length(vc)),vc),plot.na=F,mark.groups=F)+xlim(range(ncdcon$embedding[,1])) +ylim(range(ncdcon$embedding[,2]))
plot_grid(plotlist=list(p1,p2,p3),nrow=1)

Annotation of the adrenal data

neww_annot <- readRDS("~pkharchenko/m/ninib/NB/neww_annot.rds")

ap2 <- ncdcon$samples$Adr
ap2$getEmbedding(type='PCA',embeddingType='tSNE',n.cores=30)

ap2$getKnnClusters(type='PCA',method=leiden.community,r=1)

ap2$getKnnClusters(type='PCA',method=leiden.community,r=1)

ap2$plotEmbedding(type='PCA',embeddingType = 'tSNE',mark.clusters = T,cex=0.2,mark.cluster.cex = 1.5)

write out p2 app

suppressMessages(library(org.Hs.eg.db))
ids <- unlist(lapply(mget(colnames(ap2$counts),org.Hs.egALIAS2EG,ifnotfound=NA),function(x) x[1]))
rids <- names(ids); names(rids) <- ids;
# list all the ids per GO category                                                            
go.env <- list2env(eapply(org.Hs.egGO2ALLEGS,function(x) as.character(na.omit(rids[x]))))
ap2$testPathwayOverdispersion(go.env,verbose=T,correlation.distance.threshold=0.95,recalculate.pca=F,top.aspects=15)

library(GO.db)
termDescriptions <- Term(GOTERM[names(go.env)]); # saves a good minute or so compared to individual lookups                                                                         
sn <- function(x) { names(x) <- x; x}
geneSets <- lapply(sn(names(go.env)),function(x) {
  list(properties=list(locked=T,genesetname=x,shortdescription=as.character(termDescriptions[x])),genes=c(go.env[[x]]))
})
adr.p2app <- make.p2.app(ap2, dendrogramCellGroups = ap2$clusters$PCA$multilevel, geneSets = geneSets,innerOrder='odPCA');
adr.p2app$serializeToStaticFast(binary.filename = 'adr.p2.app.bin',verbose=T)

Get markers

de <- ap2$getDifferentialGenes(type='PCA',n.cores=30,append.auc=TRUE,z.threshold=0,upregulated.only=T)

source("~/m/p2/conos/R/plot.R")
fac <- ap2$clusters$PCA[[1]]
pp <- plotDEheatmap(ap2,fac,de,n.genes.per.cluster = 10 ,show.gene.clusters=T,row.label.font.size = 6)
pdf(file='adr.heatmap.pdf',width=7,height=15); print(pp); dev.off();
pp

ncdconA <- readRDS("~pkharchenko/m/ninib/NB/ncdconA.rds")

cf <- unlist(list(as.factor(setNames(rep('NB',length(nfac)),names(nfac))),neww_annot))
p1 <- ncdconA$plotGraph(alpha=0.05,size=0.8,groups=cf,plot.na=F,font.size=3,palette=function(n) c('gray95',rainbow(n-1,v=1)))
fac <- ap2$clusters$PCA[[1]]
cf <- unlist(list(as.factor(setNames(rep('NB',length(nfac)),names(nfac))),fac))
p3 <- ncdconA$plotGraph(alpha=0.05,size=0.8,groups=cf,plot.na=F,font.size=3,palette=function(n) c('gray95',rainbow(n-1,v=1)))
p2 <- ncdconA$plotGraph(alpha=0.03,size=0.5,groups=nfac,plot.na=F,font.size=3)
plot_grid(plotlist=list(p1,p3,p2),nrow=1)

Linear deviation tests

p1 <- ncdcon$plotGraph(alpha=0.03,size=0.5,groups=ncdcon$getDatasetPerCell(),plot.na=F,mark.groups=F,palette=sample.pal)
p2 <- ncdcon$plotGraph(alpha=0.03,size=0.5,groups=annot,plot.na=F,mark.groups=T,palette=annot.palf)
#p3 <- ncdcon$plotGraph(alpha=0.03,size=0.5,plot.na=F,mark.groups=T,clustering='leiden')
p3 <- ncdcon$plotGraph(alpha=0.03,size=0.5,groups=dannot,plot.na=F,mark.groups=T,palette=dannot.palf)
plot_grid(plotlist=list(p1,p2,p3),nrow=1)

df <- data.frame(ncdcon$embedding)
df$type <- annot[rownames(df)]
colnames(df) <- c('x','y','t')
ggplot(df,aes(x,y,color=t))+geom_point(size=0.1,alpha=0.1) + geom_vline(xintercept = c(-18,-13,-3,10),linetype=3)

cells <- intersect(names(annot)[annot %in% c("Mesenchymal",'SCP-like','Noradrenergic')], rownames(ncdcon$embedding)[ncdcon$embedding[,1] <10 & ncdcon$embedding[,1] > -13 & ncdcon$embedding[,2]< -17])

mes.cells <- intersect(names(annot)[annot %in% c("Mesenchymal",'SCP-like','Noradrenergic')], rownames(ncdcon$embedding)[ncdcon$embedding[,1] < -3 & ncdcon$embedding[,1] > -18 & ncdcon$embedding[,2]< -17]) 

Principal curve fit

fit.pc.pseudotime <- function(con, cells,ndims=10,embedding=NULL, return.details=F) {
  require(princurve)
  if(is.null(embedding)) {
    if(is.null(con$misc$embeddings) || is.null(con$misc$embeddings$linfit)) {
      old.emb <- con$embedding;
      embedding <- con$misc$embeddings$linfit <- con$embedGraph(method='largeVis',target.dims=ndims)
      con$embedding <- old.emb;
    } else {
      embedding <- con$misc$embeddings$linfit
    } 
  }
  embedding <- embedding[rownames(embedding) %in% cells,]
  x <- principal_curve(embedding)
  if(return.details) { 
    # identify closest cell for each principal point
    en <- conos:::n2CrossKnn(x$s,embedding,3,verbose=F,indexType='L2')
    colnames(en) <- names(x$lambda); rownames(en) <- rownames(embedding)
    en@x <- exp(mean(en@x)/en@x/10)
    en <- t(en)/colSums(en)
    x$en <- en;
    x$embedding <- embedding;
    return (x)
  }
  x$lambda
}

set.seed(0)
pcurve <- fit.pc.pseudotime(ncdcon,cells,ndims=5,return.details=T)
pt <- pcurve$lambda

pt.pos <- ncdcon$embedding[names(pt)[order(pt)],] %>% zoo::rollapply(500,median) %>% data.frame()
colnames(pt.pos) <- c('x','y');
pt.pos$sd <- ncdcon$embedding[names(pt)[order(pt)],] %>% zoo::rollapply(500,sd) %>% apply(1,sum)

Aspect ration difference between the selected subset and the actual one

Reduce(eval('/'),apply(ncdcon$embedding,2,function(x) diff(range(x))))/Reduce(eval('/'),apply(ncdcon$embedding[cells,],2,function(x) diff(range(x))))

# estimate principal curve positions
p1 <- ncdcon$plotGraph(alpha=0.2,size=0.2,colors=pt-mean(pt),plot.na=F,gradient.range.quantile=0.95) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.grid.major = element_blank())+
  geom_line(data=pt.pos[1:3.2e3,],aes(x,y),alpha=0.8,size=0.5) 
p1

Focused bridge and CNV picture for fig.2

p0 <- ncdcon$plotGraph(alpha=0.2,size=0.2,groups=annot[cells],palette=annot.pal,mark.groups=F,plot.na=F,gradient.range.quantile=0.95,raster=T,raster.width=1.5,raster.height=1.5) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) + 
  #geom_label(x=-Inf,y=Inf,vjust=1,hjust=0,label='clusters',data=data.frame(x=c(1)),label.size=0) +
  theme(panel.grid.major = element_blank())

col <- cnv.sc>0.7; 
emb <- ncdcon$embedding
ncdcon$embedding <- emb[names(col)[order(col)],]
p3 <- ncdcon$plotGraph(alpha=0.05,size=0.3,groups=col[cells],mark.groups=F,raster=T,palette=c("TRUE"='red',"FALSE"='gray90'),plot.na=F,raster.height=1.5,raster.width=1.5)+
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) + 
  geom_label(x=-Inf,y=-Inf,vjust=0,hjust=0,label='(CNV)',data=data.frame(x=c(1)),label.size=0) +
  theme(panel.grid.major = element_blank())
ncdcon$embedding <- emb;
#pdf(file='cnv_overview.pdf',width=4,height=4); print(p3); dev.off();
plot_grid(plotlist=list(p0,p3),ncol=1)

gl <- rev(c("S100B","SOX10"))
pl <- lapply(gl,function(g) ncdcon$plotGraph(colors=conos:::getGeneExpression(ncdcon,g)[cells],alpha=0.3,size=0.1,plot.na=F,gradient.range.quantile=0.99,raster=T,raster.width=1.5,raster.height=1.5)+
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) + 
  geom_label(x=-Inf,y=-Inf,vjust=0,hjust=0,label=g,data=data.frame(x=c(1)),label.size=0) +
  theme(panel.grid.major = element_blank()))
pp <- plot_grid(plotlist=c(list(p0,p3),pl),ncol=1)
pp

# returns normalized mean residuals per gene within the cel3 population, relative to a simple linear model mixing cel1 and cel2 populations
# cel1,cel2,cel3 - cell name vectors
# mat - matrix of molecule counts (raw counts - not normalized)
lin.dev <- function(cel1,cel2,cel3,mat) {
  mat <- mat[rownames(mat) %in% c(cel1,cel2,cel3),,drop=F]
  cf <- setNames(as.factor(rep(c('g1','g2'),c(length(cel1),length(cel2)))),c(cel1,cel2))
  # calculate aggregate count profiles
  cm <- conos:::collapseCellsByType(mat,cf)
  cm1 <- cm[1,]; cm2 <- cm[2,]
  # fit linear model on the bridge
  x <- as.matrix(t(mat[cel3,]))
  y <- lm(x~cm1+cm2-1)
  # report residuals
  #rs <- rstudent(y)
  # calculate pearson residual using sd from all three populations (instead of a standard studentized resiudal)
  rs <- residuals(y,'response')
  rss <- rowMeans(rs)
  rss <- rss/apply(mat,2,sd) # all three populations
  rss <- rss[is.finite(rss)]
  rss <- rss[order(rss,decreasing=T)]
}

# a convenience wrapper, starting with a pseudotime, conos object and the start/end region pseudotime ranges (reg1,reg2 are tuples of time values)
lin.dev.pt <- function(pt,con,reg1,reg2) {
  # determine cells associated with the two regions
  cel1 <- names(pt)[pt>=reg1[1] & pt<=reg1[2]]
  cel2 <- names(pt)[pt>=reg2[1] & pt<=reg2[2]]
  cel3 <- names(pt)[pt>reg1[2] & pt<reg2[1]]
  
  # construct a joint count matrix
  mat <- conos:::rawMatricesWithCommonGenes(con)
  mat <- do.call(rbind,lapply(mat,function(x) x[rownames(x) %in% c(cel1,cel2,cel3),,drop=F]))
  lin.dev(cel1,cel2,cel3,mat)
}

For illustration

For the noradrenergic side: PHOX2B, LMO1, TH For the mesenchymal side: PRRX1

rss3 <- lin.dev.pt(ptr,ncdcon,c(1.5e3,2.0e3),c(2.6e3,6e3))

gl <- rev(c("PHOX2B","PRRX1","SOX10"))
pl <- lapply(gl,function(g) ncdcon$plotGraph(colors=conos:::getGeneExpression(ncdcon,g)[cells],alpha=0.3,size=0.1,plot.na=F,gradient.range.quantile=0.99,raster=T,raster.width=1.5,raster.height=1.5)+
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) + 
  geom_label(x=-Inf,y=-Inf,vjust=0,hjust=0,label=g,data=data.frame(x=c(1)),label.size=0) +
  theme(panel.grid.major = element_blank()))
pp <- plot_grid(plotlist=c(list(p0),pl),ncol=1)
pp

pdf(file='nb.fork.genes2.pdf',width=1,height=1*0.94*4); print(pp); dev.off();

ncdcon$plotGraph(colors=conos:::getGeneExpression(ncdcon,'SOX10')[cells],alpha=0.3,size=0.1,plot.na=F,gradient.range.quantile=0.999,raster=T,raster.width=1.5,raster,height=1.5) +geom_label(x=-Inf,y=-Inf,vjust=0,hjust=0,label='SOX10',inherit.aes =F,data=data.frame(x=c(1)),label.size=0)

Test for associated genes

mat <- ncdcon$getJointCountMatrix()
mat <- t(mat[names(pt)[order(pt)[1:3.2e3]],])

# F test comparing glm models of expression with and without pseudotime
test.associated.genes <- function(pt,mat,spline.df=3,n.cores=32) {
  mat <- mat[,colnames(mat) %in% names(pt)]
  df <- data.frame(do.call(rbind,mclapply(setNames(1:nrow(mat),rownames(mat)),function(i) {
    exp <- mat[i,]
    sdf <- data.frame(exp=exp,t=pt[colnames(mat)])
    # model
    m <- mgcv::gam(exp~s(t,k=spline.df),data=sdf,familly=gaussian())
    # background
    m0 <- mgcv::gam(exp~1,data=sdf,familly=gaussian())
    fstat <- 0; 
    if(m$deviance>0) {
      fstat <- (deviance(m0) - deviance(m))/(df.residual(m0)-df.residual(m))/(deviance(m)/df.residual(m))
    }
    pval <-  pf(fstat,df.residual(m0)-df.residual(m),df.residual(m),lower.tail = FALSE);
    return(c(pval=pval,A=diff(range(predict(m)))))
  },mc.cores=n.cores,mc.preschedule=T)))
  df$gene <- rownames(mat)
  df$fdr <- p.adjust(df$pval);
  df
}

pt.genes <- test.associated.genes(pt,mat)
pt.genes <- arrange(pt.genes,pval)

Heatmaps

rescale.and.center <- function(x, center=F, max.quantile=0.99) {
  mx <- quantile(abs(x),max.quantile) # absolute maximum
  if(mx==0) mx<-max(abs(x)) # in case the quantile squashes all the signal
  x[x>mx] <- mx; x[x< -1*mx] <- -1*mx; # trim
  if(center) x <- x-mean(x) # center
  x/max(abs(x)); # scale
}

gns <- pt.genes$gene[pt.genes$fdr<1e-5 & pt.genes$A>1e-2]
gns <- pt.genes$gene[pt.genes$A>1e-2][1:1000]

#gns <- pt.genes$gene[pt.genes$fdr>1e-2 & pt.genes$A>1e-2]

gm <- zoo::rollapply(as.matrix(t(mat[rev(gns),])),30,mean,align='left',partial=T) %>% apply(2,rescale.and.center,max.quantile=1-1e-3) %>% t
colnames(gm) <- colnames(mat)
# cluster
#gm.hcl <- hclust(as.dist(2-cor(t(gm))),method='ward.D2')
gm.hcl <- hclust(dist(gm),method='ward.D2')
gm.hcl.clust <- cutree(gm.hcl,40)
gm <- gm[gm.hcl$order,]

require(ComplexHeatmap)
ha = HeatmapAnnotation(
  cells=annot[colnames(gm)],  
  #patient=ncdcon$getDatasetPerCell()[colnames(gm)],
  # specify colors
  col = list(cells=annot.pal,
             patient=sample.pal
             #pseudotime=circlize::colorRamp2(c(-1, 0, 1), c('darkgreen','grey90','orange'))
  ),
  border = T
)
ra <- ComplexHeatmap::HeatmapAnnotation(df=data.frame(clust=as.factor(gm.hcl.clust[rownames(gm)])),which='row',show_annotation_name=FALSE, show_legend=FALSE, border=T)

hm <- Heatmap(gm, cluster_columns = F, cluster_rows = F,show_column_names = F,border=T, top_annotation = ha, name='expression', show_heatmap_legend = F, show_row_dend = F, show_row_names=F, left_annotation = ra)

labeled.genes <- rownames(gm)[round(seq(1,nrow(gm),length.out = 70))]; 
hm + ComplexHeatmap::rowAnnotation(link = ComplexHeatmap::anno_mark(at = match(labeled.genes,rownames(gm)), labels = labeled.genes, labels_gp = grid::gpar(fontsize = 7)))

gns2 <- c('PHOX2B','RGS5','MLLT11',"STMN2", 'PCBP2', 'MDK','HMGB1','SRP14', #'CALM2',
         #'NLRP1',
         'ABCA8','FXYD1','CDH19','ERBB3','S100B','SOX10','MPZ','S100A10','B2M','IFITM3', # 'CNN3',
         'PODXL','NPNT', 'TNNT2','FGF1',
         'MAFF','KLF4','MYC','FOSB',
         'PRRX1','CALD1','LUM','PDGFRA','COL1A2','BGN')
#gns2 <- rownames(gm)[Reduce(seq,match(c("ANXA1","PSAP"),rownames(gm)))]
#gns2 <- rownames(gm)[Reduce(seq,match(c("TMSB15A","APC"),rownames(gm)))]
gm2 <- zoo::rollapply(as.matrix(t(mat[rev(gns2),])),30,mean,align='left',partial=T) %>% apply(2,rescale.and.center,max.quantile=1-5e-3) %>% t
#gm <- apply(mat[gns,],1,rescale.and.center) %>% zoo::rollapply(10,mean,align='left',partial=T) %>% t
colnames(gm2) <- colnames(mat)

require(ComplexHeatmap)
ha = HeatmapAnnotation(
  cells=annot[colnames(gm2)],  
  #patient=ncdcon$getDatasetPerCell()[colnames(gm2)],
  # specify colors
  col = list(cells=annot.pal,
             patient=sample.pal
             #pseudotime=circlize::colorRamp2(c(-1, 0, 1), c('darkgreen','grey90','orange'))
  ),
  border = T
)

hm2 <- Heatmap(gm2, cluster_columns = F, cluster_rows = F,show_column_names = F,border=T, top_annotation = ha, name='expression', show_heatmap_legend = F, show_row_dend = F, row_names_gp = grid::gpar(fontsize = 8),use_raster=T, raster_device = "CairoPNG")
hm2

pdf(file='small.bridge.heatmap.pdf',width=4,height=4); hm2; dev.off()

null device 
          1 

Cairo::CairoPNG(file='small.bridge.heatmap.png',width=400,height=400); hm2; dev.off()

null device 
          1 

Combined large heatmap

gns <- pt.genes$gene[pt.genes$A>1e-2][1:1000]
gns <- unique(c(gns,gns2))

gm <- zoo::rollapply(as.matrix(t(mat[rev(gns),])),30,mean,align='left',partial=T) %>% apply(2,rescale.and.center,max.quantile=1-1e-3) %>% t
colnames(gm) <- colnames(mat)
# cluster
#gm.hcl <- hclust(as.dist(2-cor(t(gm))),method='ward.D2')
gm.hcl <- hclust(dist(gm),method='ward.D2')
gm.hcl.clust <- cutree(gm.hcl,40)
gm <- gm[gm.hcl$order,]

require(ComplexHeatmap)
ha = HeatmapAnnotation(
  cells=annot[colnames(gm)],  
  #patient=ncdcon$getDatasetPerCell()[colnames(gm)],
  # specify colors
  col = list(cells=annot.pal,
             patient=sample.pal
             #pseudotime=circlize::colorRamp2(c(-1, 0, 1), c('darkgreen','grey90','orange'))
  ),
  border = T
)
ra <- ComplexHeatmap::HeatmapAnnotation(df=data.frame(clust=as.factor(gm.hcl.clust[rownames(gm)])),which='row',show_annotation_name=FALSE, show_legend=FALSE, border=T)

hm <- Heatmap(gm, cluster_columns = F, cluster_rows = F,show_column_names = F,border=T, top_annotation = ha, name='expression', show_heatmap_legend = F, show_row_dend = F, show_row_names=F,use_raster=T, raster_device = "CairoPNG")

labeled.genes <- rownames(gm)[round(seq(1,nrow(gm),length.out = 70))]; #
labeled.genes <- gns2;
hm <- hm + ComplexHeatmap::rowAnnotation(link = ComplexHeatmap::anno_mark(at = match(labeled.genes,rownames(gm)), labels = labeled.genes, labels_gp = grid::gpar(fontsize = 7)))
hm

pdf(file='full.bridge.heatmap.pdf',width=3.5,height=5); hm; dev.off()
CairoPNG(file='full.bridge.heatmap.png',width=350,height=500); hm; dev.off()

Pseudotime tests

ptr <- rank(pt) # work on ranks - they're more stable
df <- data.frame(pt=ptr,annot=annot[names(ptr)])
ggplot(df,aes(x=pt,color=annot)) + geom_density() + geom_vline(xintercept = c(1.2e3,1.5e3,2.0e3,2.6e3))

Test for deviations from linear interpolation between two clusters on the principal curve

# returns normalized mean residuals per gene within the cel3 population, relative to a simple linear model mixing cel1 and cel2 populations
# cel1,cel2,cel3 - cell name vectors
# mat - matrix of molecule counts (raw counts - not normalized)
lin.dev <- function(cel1,cel2,cel3,mat) {
  mat <- mat[rownames(mat) %in% c(cel1,cel2,cel3),,drop=F]
  cf <- setNames(as.factor(rep(c('g1','g2'),c(length(cel1),length(cel2)))),c(cel1,cel2))
  # calculate aggregate count profiles
  cm <- conos:::collapseCellsByType(mat,cf)
  cm1 <- cm[1,]; cm2 <- cm[2,]
  # fit linear model on the bridge
  x <- as.matrix(t(mat[cel3,]))
  y <- lm(x~cm1+cm2-1)
  # report residuals
  #rs <- rstudent(y)
  # calculate pearson residual using sd from all three populations (instead of a standard studentized resiudal)
  rs <- residuals(y,'response')
  rss <- rowMeans(rs)
  rss <- rss/apply(mat,2,sd) # all three populations
  rss <- rss[is.finite(rss)]
  rss <- rss[order(rss,decreasing=T)]
}

# a convenience wrapper, starting with a pseudotime, conos object and the start/end region pseudotime ranges (reg1,reg2 are tuples of time values)
lin.dev.pt <- function(pt,con,reg1,reg2) {
  # determine cells associated with the two regions
  cel1 <- names(pt)[pt>=reg1[1] & pt<=reg1[2]]
  cel2 <- names(pt)[pt>=reg2[1] & pt<=reg2[2]]
  cel3 <- names(pt)[pt>reg1[2] & pt<reg2[1]]
  
  # construct a joint count matrix
  mat <- conos:::rawMatricesWithCommonGenes(con)
  mat <- do.call(rbind,lapply(mat,function(x) x[rownames(x) %in% c(cel1,cel2,cel3),,drop=F]))
  lin.dev(cel1,cel2,cel3,mat)
}

rss <- lin.dev.pt(ptr,ncdcon,c(0,1.2e3),c(2.6e3,6e3))

rss2 <- lin.dev.pt(ptr,ncdcon,c(0,1.2e3),c(1.5e3,2.0e3))

rss3 <- lin.dev.pt(ptr,ncdcon,c(1.5e3,2.0e3),c(3e3,6e3))

pdf(file='mesenchymal.bridge.residuals.pdf',width=5,height=5); print(p1); dev.off();

null device 
          1 

sim2: PODXL … NPNT, TNNT2 NPHS1, NDNF, MME, GRN? MYC, ZFAND5, MAFF show nice foci close to the bridge on the fibroblast side

ncdcon$plotGraph(alpha=0.3,size=0.3,gene='PHOX2B',plot.na=F,gradient.range.quantile=0.99)

Run GSEA

source("~/keith/me3/gosim.r")
load("~/keith/me3/org.Hs.GOenvs.RData")

pc <- rss
gsl <- ls(env=org.Hs.GO2Symbol)
gsl.ng <- unlist(lapply(sn(gsl),function(go) sum(get(go,env=org.Hs.GO2Symbol) %in% names(pc))))
gsl <- gsl[gsl.ng>=10 & gsl.ng<=2e3]

val <- rss
gsea.rss.p1 <- iterative.bulk.gsea(values=val,set.list=lapply(sn(gsl),get,env=org.Hs.GO2Symbol),power=1,mc.cores=32)

using simple PPT

require(crestree)

Start with a high-dimensional embedding

old.emb <- ncdcon$embedding;
hiemb <- ncdcon$misc$embeddings$linfit <- con$embedGraph(method='largeVis',target.dims=4)
ncdcon$embedding <- old.emb;

restrict to the cells of interest

hiemb <- hiemb[rownames(hiemb)%in% cells,]
emb <- ncdcon$embedding[rownames(ncdcon$embedding) %in% cells,]

ptree <- ppt.tree(X=t(hiemb),emb=NA,lambda=100,sigma=0.5,metrics='euclidean',M=300,plot=F,output = F)

plotppt(ptree,emb,cex.tree=0.1,cex.main=0.2,lwd.tree=1,tips=T)

ptree <- setroot(ptree,259)

ERBB3 and NRG1

gns <- c("ERBB3","NRG1"); 
p0 <- ncdcon$plotGraph(alpha=0.1,size=0.1,groups=annot,palette=annot.pal,plot.na=F,raster=T,raster.width=4,raster.height=4)
pl <- lapply(gns,function(g) ncdcon$plotGraph(alpha=0.2,size=0.1,gene=g,plot.na=F,gradient.range.quantile=0.99,raster=T,raster.width=4,raster.height=4) +geom_label(x=Inf,y=Inf,vjust=1,hjust=1,label=g,data=data.frame(x=c(1)),label.size=0))
pp <- plot_grid(plotlist=c(list(p0),pl),nrow=1)
pdf("erbb3_nrg1.pdf",width=9,height=3); print(pp); dev.off();

png 
  2 

pp

Mesenchymal

Aspect ration difference between the selected subset and the actual one

Reduce(eval('/'),apply(ncdcon$embedding,2,function(x) diff(range(x))))/Reduce(eval('/'),apply(ncdcon$embedding[mes.cells,],2,function(x) diff(range(x))))

[1] 0.5236925

p0 <- ncdcon$plotGraph(alpha=0.2,size=0.2,groups=annot[mes.cells],palette=annot.pal,mark.groups=F,plot.na=F,gradient.range.quantile=0.95,raster=T,raster.width=2,raster.height=1) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) + 
  #geom_label(x=-Inf,y=Inf,vjust=1,hjust=0,label='clusters',data=data.frame(x=c(1)),label.size=0) +
  theme(panel.grid.major = element_blank())
p0

ncdcon$plotGraph(alpha=0.5,size=0.5,colors=conos:::getGeneExpression(ncdcon,"PODXL")[mes.cells],mark.groups=F,plot.na=F,gradient.range.quantile=0.99,raster=T,raster.width=1.5,raster.height=1.5)+
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) + 
  #geom_label(x=-Inf,y=Inf,vjust=1,hjust=0,label='clusters',data=data.frame(x=c(1)),label.size=0) +
  theme(panel.grid.major = element_blank())

mult <- 1.5; pdf(file='mes.expr.pdf',width=1*mult,height=length(pl)*0.58*mult); print(pp); dev.off()

null device 
          1